The following files will be used in this lab, all available on Canvas:
You will need to download these files from Canvas, and move them from your Downloads folder to the datafiles folder that you made previously. Make sure to unzip the zip files so that R can access the content. Note that on Windows, you will need to right-click on the file and select ‘Extract files’.
You will also need to make sure the following packages are installed on your computer:
pkgs <- c("mapview",
"raster",
"RColorBrewer"
"sf",
"tmap",
"viridis")
install.packages(pkgs)
library(mapview)
library(raster)
library(RColorBrewer)
library(sf)
library(tmap)
library(viridis)
tmaptmap is built on top of ggplot2 and works in a similar way by building a series of layers with map geometries and elements. We start by using tm_shape() to identify the spatial object to be used, and then geometries are added, including filled polygons, borders, legends, etc.
We’ll start by making some maps based on the NY shapefile. Start by loading this, then extract only the polygons belonging to Syracuse.
NY8 <- st_read("../datafiles/NY_data/NY8_utm18.shp")
## Reading layer `NY8_utm18' from data source
## `/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/NY_data/NY8_utm18.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## Projected CRS: WGS 84 / UTM zone 18N
Syracuse <- NY8[NY8$AREANAME == "Syracuse city", ]
First, let’s make a simple map showing the polygon outlines using tm_borders():
tm_shape(Syracuse) + tm_borders()
The function tm_fill() can be used to define the polygon fill color based on an attribute in the Syracuse data set (POP8). Note that this automatically adds a legend within the frame of the figure:
tm_shape(Syracuse) +
tm_borders() +
tm_fill("POP8")
The color scale can be changed by setting the palette argument in tm_fill(). This includes the ColorBrewer scales described above,and the different intervals. For example, to use the ‘Greens’ palette with percentile breaks:
tm_shape(Syracuse) +
tm_borders() +
tm_fill("POP8", palette = "viridis", style = "quantile")
Other map elements can be added. Here we add a longitude/latitude graticule with tm_graticules(), a north arrow with tm_compass(), and a line of text showing the date the map was made with tm_credits().
tm_shape(Syracuse) +
tm_graticules(col = "lightgray") +
tm_borders() +
tm_fill("POP8", palette = "Greens", style = "quantile") +
tm_compass(position = c("left", "bottom")) +
tm_credits("2019-10-19", position = c("right", "top"))
We’ll next make some maps with the Western North American site data. As before, we load the data, then convert to an sf object. We also load the shapefile of country outlines:
wna_climate <- read.csv("../datafiles/WNAclimate.csv")
wna_climate <- st_as_sf(wna_climate,
coords = c("LONDD", "LATDD"),
crs = 4326)
countries <- st_read("../datafiles/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp",
quiet = TRUE)
Individual symbols can be plotted on a color scale using tm_symbols.
tm_shape(wna_climate) + tm_symbols(col = "Jan_Tmp")
This takes the same arguments as tm_fill() for the color palette. We’ll use a red to blue color scale from RColorBrewer. The minus sign before the palette name reverses the order of the colors. As there is a large amount of overlap between the sites, we also add an alpha level to make the symbols transparent.
tm_shape(wna_climate) +
tm_symbols(col = "Jan_Tmp", alpha = 0.5, palette = "-RdBu")
We’ll next add country boundaries from the Natural Earth shapefile loaded earlier. Note that as this is a different spatial object, we have to use tm_shape() a second time to reference this, then use tm_borders() to add the lines.
tm_shape(wna_climate) +
tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
tm_shape(countries) +
tm_borders(col = "gray")
Note that here we’ve plotted the symbols first to set the plotting window. Unfortunately, this means that the country lines are added on top of the symbols. If we add the country outlines first, then the window is set to the globe.
tm_shape(countries) +
tm_borders(col = "gray") +
tm_shape(wna_climate) +
tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu")
To get around this, we can set the bounding box for the countries object using the bounding box of the point data. We’ll also move the legend to a better position in the map, and add a title
tm_shape(countries, bbox = st_bbox(wna_climate)) +
tm_borders(col = "gray") +
tm_shape(wna_climate) +
tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
tm_layout(main.title = "W. North America Climate",
legend.position = c("left", "bottom"))
We can also use one of tmap’s preset themes with tm_style():
tm_shape(countries, bbox = st_bbox(wna_climate)) +
tm_style("cobalt") +
tm_borders(col = "gray") +
tm_shape(wna_climate) +
tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
tm_layout(main.title = "W. North America Climate",
main.title.size = 1,
legend.position = c("left", "bottom"))
Finally, we can write this map to a file using tmap_save(). To do this, assign the output of all the tmap functions to an R object. Then specify this object and a file name. I’ve used pdf output here, but you can also use png, tif or most other common format extensions.
map1 <- tm_shape(countries, bbox = st_bbox(wna_climate)) +
tm_borders(col = "gray") +
tm_shape(wna_climate) +
tm_symbols(col = "Jan_Tmp", alpha = 0.75, palette = "-RdBu") +
tm_layout(main.title = "W. North America Climate",
legend.position = c("left", "bottom")) +
tm_style("cobalt")
tmap_save(wna_map, "wna_jan_tmp.pdf")
tmap includes a similar facetting function to ggplot2. This allows you to make several different versions of the same map split by some grouping factor. A good example of this is where you have repeated measures over time for a set of locations, and you would like to display each time periods separately. You can make all the individual maps, and then add them together with tmap_arrange(), or use tm_facets(). To illustrate this, we’ll load a built in dataset (metro) that has population sizes for major metropolitan centers
data(metro)
head(metro)
## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -64.18105 ymin: -34.60508 xmax: 69.17246 ymax: 36.7525
## Geodetic CRS: WGS 84
## name name_long iso_a3 pop1950 pop1960 pop1970 pop1980
## 2 Kabul Kabul AFG 170784 285352 471891 977824
## 8 Algiers El Djazair (Algiers) DZA 516450 871636 1281127 1621442
## 13 Luanda Luanda AGO 138413 219427 459225 771349
## 16 Buenos Aires Buenos Aires ARG 5097612 6597634 8104621 9422362
## 17 Cordoba Cordoba ARG 429249 605309 809794 1009521
## 25 Rosario Rosario ARG 554483 671349 816230 953491
## pop1990 pop2000 pop2010 pop2020 pop2030 geometry
## 2 1549320 2401109 3722320 5721697 8279607 POINT (69.17246 34.52889)
## 8 1797068 2140577 2432023 2835218 3404575 POINT (3.04197 36.7525)
## 13 1390240 2591388 4508434 6836849 10428756 POINT (13.23432 -8.83682)
## 16 10513284 12406780 14245871 15894307 16956491 POINT (-58.40037 -34.60508)
## 17 1200168 1347561 1459268 1562509 1718192 POINT (-64.18105 -31.4135)
## 25 1083819 1152387 1298073 1453814 1606993 POINT (-60.63932 -32.94682)
This is an sf object, with a set of columns with population size for each decade between 1950 and 2030. To plot this as a facetted map, we first need to convert this to long format with one column for area, one for population size and one for year. We’ll do this imply by concatenating the different populations sizes, and repeating the location names and year, then put it all into a temporary data frame:
pop <- c(metro$pop1950, metro$pop1960, metro$pop1970,
metro$pop1980, metro$pop1990, metro$pop2000)
name <- rep(metro$name, 6)
year <- rep(seq(1950, 2000, by = 10), each = nrow(metro))
tmp <- data.frame(pop, name, year)
Now we’ll merge this with the original sf object using the location names:
metro_long <- merge(metro, tmp, by = "name")
The merge() function acts in a similar way to SQL’s JOIN function. The name= argument allows you to specify which columns are to be used to merge.
Now we can go ahead and plot this. We’ll use symbol size to represent the population size and facet by year:
tm_shape(metro_long) +
tm_facets(by = "year") +
tm_symbols("pop", col = "darkred") +
tm_shape(countries) +
tm_borders() +
tm_layout(main.title = "Metropolitan center population")
In the final example, we’ll make figures using the global air temperature dataset. Start by re-reading the data (we also change the name of the layer, to reduce the amount of typing later on…).
air_temp <- rotate(raster("../datafiles/air.mon.ltm.nc", varname = "air"))
air_temp
## class : RasterLayer
## dimensions : 73, 144, 10512 (nrow, ncol, ncell)
## resolution : 2.5, 2.5 (x, y)
## extent : -178.75, 181.25, -91.25, 91.25 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995
## values : -37.77505, 33.58953 (min, max)
## z-value : 0000-12-30
names(air_temp) <- "jan_tmp"
proj4string(air_temp) <- CRS("+init=epsg:4326")
Now let’s make a plot using tm_raster(), which again takes similar options for color palettes. We’ll also add the country borders.
tm_shape(air_temp) +
tm_raster(col = "jan_tmp", style = "fisher", palette = "-RdBu") +
tm_shape(countries) +
tm_borders()
We can improve this a little by moving the color legend outside of the plotting area. We’ll increase the number of color classes to 9, and add a histogram showing the frequency of different values.
tm_shape(air_temp) +
tm_raster("jan_tmp",
style = "fisher",
palette = "-RdBu",
legend.hist = TRUE,
n = 9) +
tm_shape(countries) +
tm_borders() +
tm_layout(legend.outside = TRUE, legend.outside.position = "left")
For helpful tips and tricks, try:
tmap_tip()
## The filter argument of tm_shape can be used to filter spatial units.
## New since tmap 2.0
##
## data(World)
## tm_shape(World, filter = World$continent == "Africa") +
## tm_polygons("HPI")
mapviewThe mapview package provides functions to produce interactive maps in R by building on the R leaflet package, which provides a low-level interface to the javascript library of the same name. mapview supports multiple spatial data types, including points, lines, polygons, and rasters.
mapview(wna_climate,
col.regions = "darkred", # fill color
color = "gray", # outline color
alpha.regions = 0.2, # fill transparency
alpha = 0.3) # outline transparency
mapview(Syracuse,
col.regions = "darkgreen",
color = "black")
fn <- system.file("extdata", "kiliNDVI.tif", package = "mapview")
kilimanjaro_NDVI <- raster::stack(fn)[[1]]
mapview(kilimanjaro_NDVI,
col.regions = viridis::viridis(n = 10)) # supply custom color palette
Here we will only review a few of the many features that mapview provides.
The default basemap called when using mapview() is CartoDB.positron. To use other third-party map tiles, we use the map.types argument in mapview(). For a list of basemap services, see here. For convenience, the leaflet package also provides a named vector of these tile services supported by the leaflet plugin. The vector is called providers and you can access elements of the vector like you would any other:
leaflet::providers[[35]]
## [1] "MapBox"
leaflet::providers$CartoDB.Positron # for the minimalists out there
## [1] "CartoDB.Positron"
mapview(Syracuse,
col.regions = "darkgreen",
color = "black",
map.types = c("Stamen.Watercolor", "Stamen.TonerLite"))
mapview provides for layering much in the same way that ggplot2 does.
Syracuse_centers <- st_centroid(Syracuse)
lyr1 <- mapview(Syracuse,
col.regions = "darkgreen",
color = "black")
lyr2 <- mapview(Syracuse_centers,
col.regions = "darkblue",
cex = 3) # size of point
lyr1 + lyr2
We can also map attributes of the spatial data to various aesthetics.
mapview(Syracuse,
zcol = "POP8",
col.regions = viridis::cividis(n = 10),
alpha.regions = 0.85)
# breweries in Franconia, Germany
franconia_breweries <- mapview::breweries
# set point size to number of different types of beer served at each brewery
mapview(franconia_breweries, cex = "number.of.types")
For whatever reason, mapview does not currently support interactive inset maps. However, these can still be added using leaflet and a little tinkering with the mapview object.
add_miniMap <- function(x) leaflet::addMiniMap(x@map)
imap <- mapview(franconia_breweries, cex = "number.of.types")
add_miniMap(imap)